#include <vector>
#include <stdio.h>
#include <math.h>


using namespace std;
typedef unsigned int uint;


struct Vector3
{
    Vector3() : x(0), y(0), z(0) { }
    Vector3(double x, double y, double z) : x(x), y(y), z(z) { }
    double x, y, z;
};
Vector3 operator+(const Vector3& a, const Vector3& b)
{
    return Vector3(a.x+b.x, a.y+b.y, a.z+b.z);
}
Vector3 operator-(const Vector3& a, const Vector3& b)
{
    return Vector3(a.x-b.x, a.y-b.y, a.z-b.z);
}
Vector3 operator*(double a, const Vector3& b)
{
    return Vector3(a*b.x, a*b.y, a*b.z);
}
Vector3 operator/(const Vector3& a, double b)
{
    return Vector3(a.x/b, a.y/b, a.z/b);
}
bool operator==(const Vector3& a, const Vector3& b)
{
    return (a.x==b.x) && (a.y==b.y) && (a.z==b.z);
}
double Length(const Vector3& a)
{
    return sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
}
Vector3 Cross(const Vector3& a, const Vector3& b)
{
    double x = a.y*b.z - a.z*b.y;
    double y = a.z*b.x - a.x*b.z;
    double z = a.x*b.y - a.y*b.x;
    return Vector3(x, y, z);
}
double Dot(const Vector3& a, const Vector3& b)
{
    return a.x*b.x + a.y*b.y + a.z*b.z;
}
Vector3 Normalize(const Vector3& a)
{
    double length = Length(a);
    return Vector3(a.x/length, a.y/length, a.z/length);
}


void Print(const char* path, const vector<Vector3>& vertices, const vector<size_t>& indices)
{
    FILE* file;
	fopen_s(&file, path, "w");

    vector<size_t> wireframeIndices;
    for (size_t i=0; i!=indices.size(); i+=3)
    {
        size_t a = indices[i+0];
        size_t b = indices[i+1];
        size_t c = indices[i+2];
        if (a < b)
        {
            wireframeIndices.push_back(a); wireframeIndices.push_back(b);
        }
        if (b < c)
        {
            wireframeIndices.push_back(b); wireframeIndices.push_back(c);
        }
        if (c < a)
        {
            wireframeIndices.push_back(c); wireframeIndices.push_back(a);
        }
    }

	double minDistance = 1;
    for (size_t i=0; i!=indices.size(); i+=3)
    {
		Vector3 a = vertices[indices[i+0]];
		Vector3 b = vertices[indices[i+1]];
		Vector3 c = vertices[indices[i+2]];
		Vector3 barycenter = (a + b + c) / 3;
		Vector3 n = Normalize(Cross(b-a, c-b));
	    double distance = Dot(n, barycenter);
		if (distance < minDistance)
			minDistance = distance;
	}
    double epsilon = 0.001f;

    size_t vertexCount = vertices.size();
    fprintf(file, "float unitSphereVertices[%u] =\n", vertexCount*3);
    fprintf(file, "{");
    for (size_t i=0; i!=vertexCount; ++i)
    {
        Vector3 v = vertices[i];
        if (i % 2 == 0)
            fprintf(file, "\n    ");
        fprintf(file, "%+ff, %+ff, %+ff,  ", v.x, v.y, v.z);
    }
    fprintf(file, "\n};\n\n");

    fprintf(file, "float encloseUnitSphereVertices[%u] =\n", vertexCount*3);
    fprintf(file, "{");
    for (size_t i=0; i!=vertexCount; ++i)
    {
        Vector3 v = (1+epsilon) * vertices[i] / minDistance;
        if (i % 2 == 0)
            fprintf(file, "\n    ");
        fprintf(file, "%+ff, %+ff, %+ff,  ", v.x, v.y, v.z);
    }
    fprintf(file, "\n};\n\n");

    size_t indexCount = indices.size();
    fprintf(file, "uint32 triangleIndices[%u] =\n", indexCount);
    fprintf(file, "{");
    for (size_t i=0; i!=indexCount; ++i)
    {
        if (i % 9 == 0)
            fprintf(file, "\n    ");
        fprintf(file, "%u, ", indices[i]);
    }
    fprintf(file, "\n};\n\n");

    size_t wireframeIndexCount = wireframeIndices.size();
    fprintf(file, "uint32 lineIndices[%u] =\n", wireframeIndexCount);
    fprintf(file, "{");
    for (size_t i=0; i!=wireframeIndexCount; ++i)
    {
        if (i % 8 == 0)
            fprintf(file, "\n    ");
        fprintf(file, "%u, ", wireframeIndices[i]);
    }
    fprintf(file, "\n};\n\n");

    fclose(file);
}


void Weld(const vector<Vector3>& src, vector<Vector3>& vertices, vector<size_t>& indices)
{
    vertices.clear();
    indices.clear();
    size_t srcCount = src.size();
    for (size_t i=0; i!=srcCount; ++i)
    {
        Vector3 v = src[i];
        size_t vertexCount = vertices.size();
        size_t idx = vertexCount;
        for (size_t j=0; j!=vertexCount; ++j)
        {
            if (v == vertices[j])
            {
                idx = j;
                break;
            }
        }
        if (idx == vertexCount)
            vertices.push_back(v);
        indices.push_back(idx);
    }
}


void WeldAndPrint(const char* filename, const vector<Vector3>& vertices)
{
    vector<Vector3> weldedVertices;
    vector<size_t> weldedIndices;
    Weld(vertices, weldedVertices, weldedIndices);
    Print(filename, weldedVertices, weldedIndices);
}


void Subdivide(const vector<Vector3>& src, vector<Vector3>& vertices)
{
    vertices.clear();
    size_t srcCount = src.size();
    size_t i = 0;
    while (i != srcCount)
    {
        Vector3 a = src[i++];
        Vector3 b = src[i++];
        Vector3 c = src[i++];
        Vector3 ab = Normalize(0.5 * (a + b));
        Vector3 bc = Normalize(0.5 * (b + c));
        Vector3 ca = Normalize(0.5 * (c + a));
        vertices.push_back(a); vertices.push_back(ab); vertices.push_back(ca);
        vertices.push_back(b); vertices.push_back(bc); vertices.push_back(ab);
        vertices.push_back(c); vertices.push_back(ca); vertices.push_back(bc);
        vertices.push_back(ab); vertices.push_back(bc); vertices.push_back(ca);
    }
}


int main()
{
    const double phi = (1.0 + sqrt(5.0)) / 2.0;

    Vector3 srcVertices[12] =
    {
        Vector3(   0,    1,  phi), Vector3(   0,   -1,  phi), Vector3(   0,    1, -phi), Vector3(   0,   -1, -phi),
        Vector3(   1,  phi,    0), Vector3(  -1,  phi,    0), Vector3(   1, -phi,    0), Vector3(  -1, -phi,    0),
        Vector3( phi,    0,    1), Vector3( phi,    0,   -1), Vector3(-phi,    0,    1), Vector3(-phi,    0,   -1),
    };
    uint srcIndices[60] =
    {
        4,  5,  0,   0,  5, 10,   0, 10,  1,   8,  0,  1,   4,  0,  8,
        4,  8,  9,   9,  8,  6,   8,  1,  6,   6,  1,  7,   6,  7,  3,
        6,  3,  9,   9,  3,  2,   2,  3, 11,   2, 11,  5,   2,  5,  4,
        9,  2,  4,  11, 10,  5,   7, 10, 11,   1, 10,  7,  11,  3,  7,
    };

    vector<Vector3> vertsA, vertsB;
    for (size_t i=0; i<60; ++i)
        vertsA.push_back(Normalize(srcVertices[srcIndices[i]]));

    WeldAndPrint("geosphere0.txt", vertsA);
    Subdivide(vertsA, vertsB);
    WeldAndPrint("geosphere1.txt", vertsB);
    Subdivide(vertsB, vertsA);
    WeldAndPrint("geosphere2.txt", vertsA);
    Subdivide(vertsA, vertsB);
    WeldAndPrint("geosphere3.txt", vertsB);
    Subdivide(vertsB, vertsA);
    WeldAndPrint("geosphere4.txt", vertsA);

    return 0;
}
